Discovery of two-dimensional binary nanoparticle superlattices using global Monte Carlo optimization

Binary nanoparticle (NP) superlattices exhibit distinct collective plasmonic, magnetic, optical, and electronic properties. Here, we computationally demonstrate how fluid-fluid interfaces could be used to self-assemble binary systems of NPs into 2D superlattices when the NP species exhibit different miscibility with the fluids forming the interface. We develop a basin-hopping Monte Carlo (BHMC) algorithm tailored for interface-trapped structures to rapidly determine the ground-state configuration of NPs, allowing us to explore the repertoire of binary NP architectures formed at the interface. By varying the NP size ratio, interparticle interaction strength, and difference in NP miscibility with the two fluids, we demonstrate the assembly of an array of exquisite 2D periodic architectures, including AB-, AB2-, and AB3-type monolayer superlattices as well as AB-, AB2-, A3B5-, and A4B6-type bilayer superlattices. Our results suggest that the interfacial assembly approach could be a versatile platform for fabricating 2D colloidal superlattices with tunable structure and properties.

NP superlattices have attracted considerable scientific interest because of their potential applications in plasmonic 1-3 , optoelectronic 4 , photonic 5 , and phononic 6,7 devices as well as in catalysis 8,9 and energy conversion 10,11 . Compared to single-component NP superlattices, binary NP superlattices (BNSLs)-composed of two types of NPs differing in size or composition-possess significantly more structural and chemical diversity. As a result, BNSLs can exhibit collective properties distinct from those of individual NPs, disordered NP structures, and single-component NP superlattices [12][13][14] . For instance, BNSLs assembled from plasmonic and non-plasmonic NPs exhibit a collective plasmonic response that can be tuned across the entire visible spectrum by varying NP size, composition and lattice symmetry 15 . Among them, 2D BNSLs are especially attractive because they could also be used as templates in colloidal lithography for fabricating nanostructured materials 16 and applications such as electroluminescence prefer the NPs to be assembled into monolayers rather than 3D structures 17,18 .
Interfaces formed between mutually immiscible fluids provide a powerful platform for assembling NPs into 2D higher-order structures 19 . However, spherical NPs tend to form hexagonally packed monolayers at fluid-fluid interfaces to maximize interparticle contacts 20 . Recently, we proposed a strategy for assembling more unusual structures from spherical NPs trapped at interfaces 21 . The approach takes advantage of the phenomenon where particles get increasingly displaced from the interfacial plane with increasing difference in the surface energy of NPs with the two fluids. Such differences could be engineered by grafting NP surfaces with ligands that have different miscibility with the two fluids. By co-assembling multiple species of grafted NPs, where species get trapped in distinct planes parallel to the interfacial plane based on how differently their grafts interact with the two fluids, we demonstrated via coarse-grained molecular dynamics (CGMD) simulations that such NPs could be assembled into distinct clusters, strings, and sheets. Among layered structures relevant to 2D BNSLs, we successfully assembled a ridged monolayer and a square-packed interdigitated bilayer. While these examples illustrate the potential of our interfacial assembly approach for fabricating 2D BSNLs, only a handful of assembly scenarios could be investigated due to the high computational cost of MD simulations, especially for superlattices with large, complex unit cells. Therefore, the full potential of this approach for self-assembling 2D BNSLs could not be explored. Furthermore, there is no guarantee in MD simulations that the final, assembled NP configuration represents the globally stable state, as configurations often get trapped in metastable states. In fact, the number of such locally stable states increases exponentially with system size, and so do the number of transition states 22,23 . The task of obtaining globally stable NP structures is thus a global optimization problem. Many optimization algorithms have been developed for molecular and crystal structure prediction [24][25][26][27][28][29][30] , but none are applicable to studying interfacial assembly of NPs. Similarly, although some computational efforts have been devoted to exploration of NP superlattices [31][32][33][34][35][36][37] , none of them have focused on 2D BNSLs formed at interfaces.
In this work, we develop an efficient global optimization approach tailored for obtaining the ground-state structure of interface-trapped NP assemblies and explore the full range of structural phases that binary systems of NPs form at fluid-fluid interfaces (Fig. 1). By examining a range of NP size ratios, interparticle interaction strengths, and surface tension differences of NPs with the liquid bilayer, we uncover many interesting monolayer, bilayer, and globular structures. By analyzing large-scale configurations of the obtained monolayer and bilayer clusters, many unusual 2D BNSLs structures are revealed. The optimization approach designed here, the assembly strategy explored here, the understandings gained here, and the structures discovered here, all represent major advancements in the synthesis and discovery of 2D materials.

Model development and structure optimization
The system of interest comprises of two species of surfacefunctionalized spherical NPs (labelled 1 and 2) trapped at a planar interface (located at position z = 0) formed between two mutually immiscible fluids (denoted by 1' for fluid layer at z < 0 and 2' for layer at z > 0) in which the two types of NPs interact favorably with opposite layers of the fluid bilayer, that is, NP 1 with fluid 1' and NP 2 with fluid 2' (Fig. 1a) c Optimization outcome for a binary 6:6 system of NPs (φ = 1, ϵ = 0:08, χ = 0:6, ω = 1) that yields a square-packed bilayer as the global minimum. The total energy (normalized by πR 2 1 γ) vs. MC steps and the various structures obtained during the optimization process are shown. The large initial drop in energy is attributed to a transition from a loosely dispersed state of NPs to a relatively dense planar structure. The structure corrects itself in subsequent steps of optimization, transitioning through increasingly less-defective square-packed bilayer structures. At the penultimate step, the algorithm obtains a square-packed bilayer with parallelaligned NP layers, before yielding the globally stable square-packed bilayer with perpendicularly aligned layers.
interparticle interactions and interfacial interactions of the two NPs with the two fluid layers. As detailed in Methods, our model treats the two types of NPs as spheres of radii R 1 and R 2 that interact with each other via short-ranged potentials U ij (ij = 11, 22, and 12 denotes interactions between like and unlike pairs of NPs) whose interaction strengths ε ij were assumed to be proportional to NP radii, consistent with the size scaling of van der Waals interactions between spherical colloids (Supplementary Fig. 1) 38 . Throughout this work, we fix R 1 = 3σ (σ represents an intrinsic length scale of the system) and vary the radius R 2 of the other NP type. To describe the interaction of NPs with the fluid bilayer, we used an analytical model we recently developed 21 for the position (z)-dependent free energy 4F i of a surface-functionalized NP. The model accounts for the surface tension γ between the fluids, which attempts to trap the NPs symmetrically about the interfacial plane so that they occlude the largest possible interfacial area, and difference in the surface energy 4γ i of the NPs with the two fluids, which attempts to expel NPs from their less favored fluid layer ( Supplementary Fig. 2). The interplay between these two effects determines the equilibrium position z m of the NPs (in the absence of interparticle interactions): where the sign in Eq. (1) is "À" for NP 1 and " + " for NP 2, and 4F i z m À Á is the interfacial free energy of the NP at this position. The total energy of a binary system of N NPs (N 1 of type 1 and N 2 of type 2, so N = N 1 + N 2 ) is then represented by the sum of NðN À 1Þ=2 pairwise interactions between the NPs and N interfacial energy contributions from individual particles.
Dimensional analysis over the energy expressions leads to four dimensionless parameters that entirely describe the assembly thermodynamics of the NPs: where φ is the size ratio of the two species of NPs, ϵ is the relative strength of interparticle interaction ε 11 to surface tension, χ is the ratio of the difference in the surface energy of NP 1 with the two fluids to the surface tension between the two fluids (this parameter also determines the equilibrium position of NPs as seen from Eq. (1)), and ω is the ratio of the differences in the surface energy of the two types of NPs with the two layers. Table 1 lists the parameter values whose combinations are explored in this study.
To locate the global minimum-energy (ground-state) configuration of NPs for each parameter combination, we used the BHMC global optimization algorithm [39][40][41] . However, the NP displacement "moves" implemented in standard BHMC algorithms are not particularly efficient for sampling NP configurations of the quasi-2D morphologies explored in this work 13,19,21 . To this end, we designed several new MC moves for improving the sampling efficiency of the BHMC algorithm ( Fig. 1b). To demonstrate the efficiency of our tailored BHMC algorithm, we present in Fig. 1c the evolution in total energy with MC steps during optimization of a representative "6:6" NP system (N 1 = N 2 = 6) whose globally stable structure is the square-packed bilayer obtained earlier through CGMD simulations of polymer-grafted NPs in a polymer bilayer 21 . The entire optimization required only 33 MC steps.

Structural phases obtained from equal-sized NPs
We first explored the assembly structures formed by a 6 : 6 system of "anti-symmetric" NPs, where the two species of NPs have the same size (φ = 1) and the same extent of preference for the opposite fluid layers (ω = 1). Figure 2a presents a "phase diagram" of globally stable structures obtained by our optimization approach for different combinations of interparticle interaction strengths ϵ and relative miscibility values χ. When both NP species interact identically with the two fluids (χ = 0), the problem reduces to an optimization of a single species of NPs. As expected, such NPs reside symmetrically at the interface and form hexagonally packed monolayers when interfacial forces dominate interparticle interactions (ϵ < 0:32). This arrangement allows the NPs to occlude the maximum possible area of the interface and gain as many interparticle interactions as possible by a planar array of particles. As χ becomes larger, the hexagonal monolayer becomes more corrugated (structures enclosed by blue boxes) because the two sets of NPs prefer to stay slightly above and below the interfacial plane. The observed ridged pattern of upward-and downward-displaced NPs allows them to minimize the separation distance between nonneighboring NPs so that they can gain more energy from long-range interparticle interactions.
Bilayers composed of two interdigitated layers of ordered NPs begin to appear upon further increase in χ in this regime dominated by interfacial forces. These bilayer structures, shown in red boxes, arise due to NPs attempting to gain more energy through their interactions with NPs of the opposite layer while trying to retain as much as possible their preferred interfacial position. A hexagonally ordered bilayer forms at χ = 0:8 in which the two layers of NPs are located at z h i = ±2:50σ relative to the interfacial plane, where Á h i denotes ensemble average over all NPs in each layer. From Eqs. (1), (5), and (6), the preferred interfacial position of such NPs without interparticle interactions is z m = ± χR = ± 2:4σ (NP radii R = 3σ), close to their eventual position in the assembled structure. Thus, the NPs sacrifice little interfacial energy in forming the hexagonal bilayer. However, at χ = 0:6, for which z m = ± 1:8σ, the NPs prefer to instead form a squareordered bilayer that incurs a lesser interfacial energy penalty, as this structure requires the NPs to be less displaced ( z h i = ± 2:16σ) from their preferred positions ( ± 1:8σ) compared to the hexagonal bilayer ( z h i = ± 2:50σ). In both cases (χ = 0:6 and 0:8), the incurred losses in interfacial free energy are easily compensated by the interparticle interaction energy gained through the formation of bilayers.
At χ = 1, where z m = ± R and 4F z m À Á = 0 (from Eq. 2), the two species of NPs would prefer to fully submerge in the fluid layers they are most compatible with and form individual close-packed globular clusters in their respective layers. However, the NPs could gain additional energy through interactions with NPs in the opposite layer. Indeed, our optimization finds that for ϵ < 0:64 the NPs form quasiglobular clusters in their respective fluid layers (shown in orange boxes), where the clusters present large facets in the direction of the interface to maximize NP-NP contacts between clusters across the interface. Unlike the bilayers discussed above, the interacting layers of NPs avoid interdigitation to avoid highly unfavorable interactions with the opposite fluid.
These results show that when interparticle interactions are weak (ϵ < 0:64), χ is the primary governing factor underlying all structural transitions. This may also be gleaned from the map of ground-state energies ( Supplementary Fig. 3), which displays much stronger variation with χ than ϵ at small values of ϵ. The map also reveals that ϵ becomes the dominating factor when ϵ>0:64. Indeed, when interparticle interactions dominate interfacial interactions, the NPs form close-packed globular clusters (shown in purple boxes), resembling the truncated triangular dipyramid recognized to be the ground-state configuration of a 12-colloid cluster in the bulk 21 . Interestingly, as χ varies, the clusters retain their truncated triangular dipyramid configuration but adjust their position and orientation with respect to the interfacial plane so that the NPs can maximize their interactions with the interface without losing interparticle interactions. Sandwiched in between monolayers and globular clusters at small χ values are hybrid morphologies (shown enclosed by cyan boxes) that arise due to compromise between interfacial and interparticle interactions, so these structures exhibit partly planar and partly globular features.
The system size used so far for exploration is not large enough to establish the true periodicity of some of the observed phases due to edge or finite size effects. For example, the square bilayer obtained from optimization appears with either parallel (ϵ = 0:04, χ = 0:6) or crossed (ϵ = 0:08, χ = 0:6) layers of NPs (Fig. 2a); at small ϵ, parallel layers are favored because they can bend at their edges to occlude more interfacial area, and at large ϵ, crossed layers are preferred as they can mediate one extra NP-NP contact compared to their parallel counterparts. Despite these differences, both structures are expected to converge to a perfectly ordered square bilayer in the limit N ! 1.
To better reveal particle arrangements expected in macroscopic versions of the predicted monolayer and bilayer structures, we applied our optimization strategy to larger systems of NPs ( Supplementary  Fig. 4), which reduces edge effects and allows periodicity to appear naturally, without imposition of periodic boundaries. Figure 2b-e presents the key structures we achieved from this analysis. We find that the ambiguous conjoined quasi-globular clusters predicted at χ = 1 now form more defined hexagonally close packed clusters in the two layers that continue to interact with each other across the interface via their large facets without NP interdigitation ( Fig. 2b and Supplementary Fig. 4a). The hexagonal bilayer structure obtained earlier retains its morphology ( Fig. 2c and Supplementary Fig. 4b), while the conditions that led to the two variants of square bilayer structures reported earlier now yield the same square-ordered bilayer ( Fig. 2d and Supplementary Fig. 4c). Lastly, the punctuated ridges observed in the hexagonal monolayer structures obtained at small χ > 0 take form and lead to zigzag patterns, arguably to maximize the long-ranged attraction between NPs as indicated earlier (Fig. 2e). In the limit of N ! 1, we expect that the planar and ridged monolayer structures would continue to grow into 2D AB-type hexagonally arranged superlattices, the bilayer structures would form 2D AB-type superlattices with square or hexagonal arrangements, and the quasi-globular clusters and the truncated triangular dipyramid would both evolve into 3D fcc or hcp superlattices.
The observed transitions in assembly morphology can be described using a simple phenomenological model (see Methods for model derivation). In this model, the total free energy (normalized by πR 2 1 γ) of a NP structure is given by the sum of its interfacial and interparticle interaction energies and approximated using Àm 1 À χ ð Þ 2 À nϵ, where m is the effective number of NPs in the structure that maintain their equilibrium positions (z m ) and n is the total number of NP-NP contacts exhibited by the structure. By determining the values of m and n for the obtained globally stable clusters (Fig. 3a-c), one can determine the χand ϵ-dependent free energies of the monolayer, bilayer, and globular phases, which we denote by F m , F b , and F g . The coexistence boundaries between the three phases can then be determined from phase equilibria conditions F m ðχ, ϵÞ = F b ðχ, ϵÞ, F m ðχ, ϵÞ = F g ðχ, ϵÞ, and F b ðχ, ϵÞ = F g ðχ, ϵÞ (Supplementary Table 1, Supplementary Fig. 5). The obtained coexistence boundaries (Fig. 3d) show excellent agreement with the earlier phase diagram. The boundaries interestingly point to the existence of a "triple point" where all three phases are in coexistence, suggesting that it may be possible to interconvert NP assemblies between the three phases through slight changes in the parameters χ and ϵ about this point.
In addition to antisymmetric NPs (ω = 1), we also examined the structural phases formed by "asymmetric" NPs, where the two NP species 1 and 2 are still similar in size but exhibit different extents of preference for opposite fluid layers (ω ≠ 1). The phase diagrams of globally stable configurations of equal-size NPs at ω = 0, ω = 0:5, and ω = 2 are presented in Supplementary Figs. 6-8. At ω = 0, NP 2 interacts identically with both fluid layers and thus would prefer to stay symmetrically at the interface. At ω = 0:5 or ω = 2, NP 2 would exhibit a preferred displacement (z m ) that is half of or twice than that of NP 1. We find that most of the structures formed by asymmetric NPs are already revealed in the phase diagram of antisymmetric NPs (Fig. 2a). The reason is because the NP arrangement in the monolayer and bilayer phases is primarily determined by the relative displacement (in the z direction) of the two types of NPs, which is given by χð1 + ωÞR (from Eqs. (1), (5), and (6)) and has been explored via changes in χ, and the structures of the globular phases are mainly determined by the interaction strength ϵ.

Structural phases obtained from different-sized NPs
Motivated by the uniqueness of some of the layered phases obtained with equal-sized NPs, we explored structures formed by unequal-sized NPs. We examined a 6 : 6 binary system of NPs and considered different NP size ratios (φ = 5=6, 2=3, 1=2, and 1=3) holding ω fixed at a value of 1 given that this parameter only weakly impacted NP phase behavior. The phase diagrams corresponding to the four size ratios are shown in Supplementary Figs. 9-12, revealing a rich variety of assembly architectures. Like equal-sized NPs, the unequal-sized NPs also form monolayers at small ϵ and χ, bilayers at small ϵ but large χ, and globular structures at large ϵ. While the globular structures exhibit compact close-packed morphologies expected to grow into 3D binary crystals akin to those formed in the bulk 42 , the monolayers and bilayers are distinct from those found in bulk solution and from those formed with equal-sized NPs. Therefore, we focused on exploring larger-scale configurations of these monolayer and bilayer phases (superlattices), which display unusual particle arrangements that would be particularly relevant to plasmonic, magnetic, optical, and electronic applications. The obtained configurations were found to be quite robust to deviations in the number ratio of the two NP species (used during structure optimization) from the particle stoichiometry stipulated by superlattice periodicity (Supplementary Fig. 13).
We begin by discussing the monolayers formed at small χ values by binary NP systems of slightly different sizes (φ = 5=6). Despite the size difference, we find that the NPs organized themselves into rough hexagonal order to maximize interparticle interactions (Supplementary Fig. 9). At χ = 0:4 and ϵ = 0:08, however, the NPs form a perfectly hexagonally ordered monolayer in which small NPs form a Kagome sublattice whose voids are filled with large NPs that form a hexagonal sublattice (Fig. 4a). At this χ value, both NP types prefer to be displaced away from the interfacial plane, with an ideal (equilibrium) spacing of χ R 1 + R 2 À Á = 0:4 × 3σ + 2:5σ ð Þ = 2:2σ between the two sets of NPs. For the small NPs to form perfect hexagons around the large NPs, geometric arguments suggest that the two types of NPs need to be displaced by a distance of 2:29σ ( Supplementary Fig. 14a). As this value is very close to the ideal spacing between the two NP types, the particles sacrifice little interfacial energy to attain the observed close-packed hexagonal arrangement of small and large NPs, which allows the NPs to gain substantial interparticle interaction energy. A larger periodic lattice of such NPs is schematically shown in Fig. 4a to illustrate how these monolayer clusters are expected to continue to grow into a 2D AB 3type NP superlattice.
Reducing the size ratio to φ = 2=3 causes the NPs to form stretched versions of hexagonally ordered monolayers (which become especially ordered at χ = 0:2), where the large NPs form a strained honeycomb lattice around dimers of small NPs (Fig. 4b).
Given that the 1:1σ spacing observed between the two NP types in this stretched monolayer is also very close to their ideal spacing of χ R 1 + R 2 À Á = 0:2 × 3σ + 2σ ð Þ= 1:0σ ( Supplementary Fig. 14b) means that this structure, like the monolayer of Fig. 4a, allows the NPs to gain more interparticle interactions without sacrificing much interfacial energy. This quasi-hexagonal arrangement of NPs is expected to grow into an AB-type 2D superlattice. Another AB-type monolayer is formed at φ = 1=2 in which the large and small NPs form intertwined square sublattices (Fig. 4c). The spacing observed between the two types of NPs in this square arrangement is 1:5σ ( Supplementary Fig. 14c), which is slightly less than the ideal spacing of χ R 1 + R 2 À Á = 0:4 × 3σ + 1:5σ ð Þ= 1:8σ. Therefore, NPs sacrifice a small amount of interfacial energy to gain the higher number of interparticle interactions afforded by the square arrangement.
Two distinct varieties of AB 2 -type monolayers form at φ = 1=3 and χ = 0:4. The first is a triangular-ordered monolayer formed by relatively strongly interacting NPs (ϵ = 0:08), where the large NPs form a hexagonally close-packed monolayer sublattice and the small NPs occupy the interstices formed between the large NPs (Fig. 4d). The other is a rhombic-ordered monolayer formed by weakly interacting NPs (ϵ = 0:02), where the large NPs arrange into rhombuses, and the small NPs form dimers and sit at the interstices of this rhombic sublattice (Fig. 4e). Based on their unit cells (marked by black lines in Fig. 4d, e), the two monolayers appear somewhat similar, as both have an underlying rhombic unit cell. However, a key difference is that the triangular monolayer shows direct contact between large NPs and no contact between small NPs, while the rhombic monolayer shows no contact between large NPs and direct contact between small NPs. Since the interaction strength between NPs is proportional to particle size, the unit cell of a triangular lattice displays more favorable interparticle interaction energy than that of a rhombic lattice. On the other hand, the triangular monolayer displays less favorable interfacial energy than its rhombic counterpart, as the small NPs exhibit more deviation from their ideal position ( z h i≈0:8σ compared to z m = χR 2 = 0:4σ) in the triangular monolayer (Supplementary Fig. 14d) than in the rhombic monolayer, where the small NPs are located almost at the ideal position ( Supplementary Fig. 14e); in both monolayers, the large NPs reside close to their ideal position, as they sacrifice much larger interfacial energy compared to small NPs due to their much larger cross-sectional area. Thus, the triangular monolayer with more favorable interparticle interactions and less favorable interfacial interactions is preferred at large values of ϵ, while the rhombic monolayer is preferred at small ϵ; a simple calculation of energy differences predicts a cut off value of ϵ ∼ 0:053 for this transition (see Methods).
Our optimization also revealed another globally stable rhombic configuration of NPs (Fig. 4f) with the same energy as the structure in Fig. 4e. In both monolayers, the large NPs form rhombuses while the small NPs sit at their interstices as dimers. The difference is that in Fig. 4e, all small-NP dimers align parallel to each other, while in Fig. 4f, the dimers in adjacent layers align at an angle of~90 . Thus, the structures could continue to grow into larger scale AAAA or ABAB layered arrangements ( Supplementary Fig. 4h), analogous to the hcp and fcc lattices formed in 3D that have the same packing fraction but different stacking: ABABAB for hcp and ABCABC for fcc. Interestingly, the NPs could also form hybrid structures with randomly mixed A and B layers, an example of which is illustrated in Fig. 4g.
Several of these predicted monolayer NP superlattices have been successfully realized in experiments using NP size ratios similar to those stipulated by our computational predictions 19,43 , as summarized in Supplementary Fig. 15. For example, the AB-type square arrangement of NPs (Fig. 4c) obtained with 2:1 size ratio was assembled at an air-liquid interface from a mixture of 28.9 nm NaFY 4 :Yb/Er + 13.4 nm Fe 3 O 4 nanocrystals (NCs) 19 (Supplementary Fig. 15c). A structure similar to our AB 2 -type triangularly ordered BNSL (Fig. 4d) predicted for a size ratio of 3:1 was assembled at an air-water interface from a mixture of 500 nm and 100 nm polystyrene NPs 44 ( Supplementary  Fig. 15d). Even the strained honeycomb (Fig. 4b) and rhombic (Fig. 4e, f) BNSLs predicted here have been experimentally realized, albeit using LaF 3 nanodisks and CdSe/CdS nanorods 43 ( Supplementary  Fig. 15g, h), which also have circular cross-sections like our spherical NPs.
When χ > 0:5, the NPs prefer to reside further away from the interfacial plane (|z m | rises with increasing χ from Eq. (1)) and their interactions with the interface become weaker (from Eq. (2), |4F z m À Á | decreases with increasing χ), causing the two sets of NPs to form bilayer configurations (Fig. 5a-d). Thus, interparticle interactions play a stronger role in dictating NP arrangement in bilayer structures than monolayers. As opposed to equal-sized NPs that form a hexagonally packed bilayer at χ = 0:8 (Fig. 2c), NPs slightly different in size (φ = 5=6) assemble into a more distinctive hybrid bilayer architecture, where the small and large NP species form layers of different symmetries (Fig. 5a). As the strength of interactions between NPs is proportional to their size, large NPs play a greater role in dictating the ground-state configuration. Consequently, the large NPs preserve their favored hexagonally ordered monolayer arrangement, while the small NPs adjust to the topography of this configuration by forming a "stretched" hexagonal monolayer with a rhombic unit cell of internal angle 106 (note that 120 angle will lead to regular hexagons). The spacing between the two layers in this hybrid bilayer is 4:6σ ( Supplementary Fig. 16a), close to the ideal value of χ R 1 + R 2 À Á = 0:8 × 3σ + 2:5σ ð Þ = 4:4σ, so the NPs sacrifice little interfacial energy to form this hybrid architecture. This structure in which small NPs sit at the interstices of the hexagonal lattice of large NPs also allows the small NPs to maximize their interactions with the large NPs while still allowing for enough contacts amongst themselves. We note however that the rhombic lattice of small NPs is inherently incompatible with the hexagonal lattice of large NPs, as we show schematically in Fig. 5a. Thus, small NPs sitting perfectly at the interstices of large NPs at one location will eventually lose this registry, leading to the formation of grain boundaries at regions of extreme incongruence ( Supplementary Figs. 16a, 4d).
Reducing the size ratio to φ = 2=3 yields another intriguing hybrid bilayer architecture (at χ = 0:8, ϵ = 0:02), where the small NPs form a pentagonal lattice (Fig. 5b) atop a hexagonal lattice of large NPs (also see Supplementary Figs. 16b, 4e). Pentagonal arrangements are rarely found in crystal structures because regular pentagons (with a 108 Random stacking of these two configurations should lead to many possible largescale hybrid arrangements; one such arrangement is shown for illustration (g). Larger-scale patterns formed by the two NPs are shown schematically at the bottom. Black lines mark the unit cells. Source data are provided as a Source Data file.
internal angle) cannot form periodic tessellations ( Supplementary  Fig. 16b). However, in our bilayer, the pentagons are slender, with one angle of 96 , two angles of 120 , and two angles of 102 . This distortion not only enables the small NPs to gain more interfacial energy (as the spacing between the two layers is close to their ideal value of χ R 1 + R 2 À Á = 0:8 × 3σ + 2σ ð Þ= 4σ; see Supplementary Fig. 16b), but it also allows the NPs to form a periodic pentagonal arrangement. Three tail-to-tail connected pentagons enclosing a triangle in the center can repeat themselves on the hexagonal lattice formed by the large NPs. Indeed, one 96 and two 102 inner angles of the three pentagons along with the 60 inner angle of the enclosed triangle add up to 360 and produce one closed node in the lattice, and three 120 inner angles of the three connected pentagons produce another kind of node. As shown in the larger scale structure in Supplementary Fig. 4e and illustrated in Fig. 5b, this arrangement is expected to continue to grow into a 2D A 3 B 5 -type superlattice.
For the same size ratio, reducing χ to 0.6 and increasing ϵ to 0.16 causes the NPs to form an equally intriguing hybrid bilayer structure comprised of a repeating pattern of hexagonal, square, and triangular motifs (Fig. 5c, also see Supplementary Figs. 16c, 4f). The large NPs still form a hexagonal lattice, but their displacement from the interfacial plane is uneven-large NPs enclosed by hexagons of small NPs stay closer to the interface than those enclosed by squares (side view, Fig. 5c). Analysis of NP spacings reveals that the small-and large-NP layers are spaced 3σ and 4:12σ apart in the hexagonal and square motifs ( Supplementary Fig. 16c). Although the hexagonal motif allows the NPs to gain the most interfacial energy (ideal spacing between the two NP types is also 3σ), the small NPs are unable to form a periodic hexagonal arrangement, as this arrangement would be incompatible with the underlying hexagonal lattice of large NPs ( Supplementary  Fig. 16c). As a result, the NPs instead adopt the multi-motif tessellation, which not only allows the NPs to gain more contacts by forming a periodic arrangement, but also maximize the appearance of hexagonal motifs in the arrangement. This arrangement leads to a 2D A 4 B 6 -type NP superlattice (Fig. 5c and Supplementary Fig. 4f).
Lastly, a periodic honeycomb lattice is formed at a size ratio of φ = 1=2 (Fig. 5d and Supplementary Fig. 4g). Here, small NPs occupy the interstices of three large NPs, culminating in an 2D AB 2 -type NP superlattice. This honeycomb lattice looks similar to the triangular lattice obtained earlier with φ = 1=3 (Fig. 4d and Supplementary   Fig. 14d), except that the small NPs in the honeycomb lattice contact each other due to their larger size, leading to closed hexagons. In addition, the spacing observed between the two NP types (2:87σ) in the honeycomb lattice is quite close to their ideal value of χ R 1 + R 2 À Á = 0:6 × 3σ + 1:5σ ð Þ= 2:7σ, implying that the small NPs can gain significant interparticle interactions without sacrificing much interfacial energy.
Among the bilayer NP superlattices reported here, the 2D AB 2type honeycomb lattice (Fig. 5d) is the only structure that has been experimentally realized, using polymer-functionalized 5 nm and 10 nm gold NPs (φ = 1=2, in agreement with our predictions) assembled at an air-water interface 45 , and 210 nm and 420 nm oppositely charged polystyrene NPs (also φ = 1=2) obtained via spin-coating method 46 ( Supplementary Fig. 15b). As far as we know, the 2D A 3 B 5 -type (Fig. 5b) and A 4 B 6 -type (Fig. 5c) superlattices predicted here have never been reported.

Discussion
This work explores the repertoire of 2D superlattices formed by a binary system of NPs at a fluid-fluid interface. Our approach of using analytical models for interparticle and interfacial interactions and an efficient Monte Carlo optimization algorithm allows us to rapidly determine the ground-state configurations of NPs trapped at an interface. By exploring NP systems of varying size ratios, interparticle interaction strengths, and miscibility differences with the two fluids, we discover a range of NP monolayer, bilayer, and globular phases. The formation of these phases is shown to arise from a competition between interfacial and interparticle interactions, and the phases boundaries are captured by an analytical model based on the effective number of NPs occupying their preferred interfacial position and the number of NP-NP contacts in the structures. By analyzing particle arrangements in large-scale configurations of these phases, we uncover a rich variety of mono-and bilayer NP superlattices. It has been demonstrated that binary NP superlattices possess plasmonic, magnetic, optical, and electronic properties distinct from those exhibited by individual NPs 18,[47][48][49] . Such collective properties depend on the symmetry and interparticle spacing of the superlattices. The diverse superlattices found in this work, combined with an understanding of their assembly revealed here, should offer opportunities for researchers to fabricate such superlattices and explore their properties and applications. Several studies have shown that NP films assembled at fluid interfaces can be successfully lifted off the interface and deposited onto a substrate 19,[50][51][52] .
While some of the predicted structures have been experimentally realized 19,43,44,46,53,54 (Supplementary Fig. 15), others such as the A 3 B 5and A 4 B 6 -type BNSLs will need to be validated and could be worthy candidates of future experimental efforts. To this end, the dimensionless parameters explored in this work would need to be linked to experimental conditions. To illustrate this process, we consider the system of silica (2) and silver (1) NPs at an air (2 0 )-water (1 0 ) interface of surface tension γ = 72:2 mJ/m 2 55 . The surface energy of silica in air (γ 22 0 ) and water (γ 21 0 ) is 287 mJ/m 2 and 340 mJ/m 2 55,56 , which yields 4γ 2 γ 21 0 À γ 22 0 = 53 mJ/m 2 . The difference in the surface energy of silver with air and water can be estimated from the contact angle of water on silver via Young's equation: 4γ 1 γ 12 0 À γ 11 0 = γ cos 64 o ≈ 31:7 mJ/m 2 57 . From Eqs. 5 and 6, we then obtain χ = 0:44 and ω = 1:67. The remaining parameter ϵ can be estimated from the measured Hamaker constant of silver in water, as given by 61.6 zJ 58 . As shown in Methods, this value yields a vdW attraction energy of ε 11 = 135 zJ between two R = 3 nm silver NPs, and then ϵ = 0:066 from Eq. (4). Therefore, according to our phase diagram, this example system would have a good chance of forming the AB-type square bilayer at size ratio φ = 1 (Fig. 2) and forming the AB 2 -type honeycomb bilayer at size ratio φ = 1=2 (Supplementary Fig. 11; Fig. 5d).
Although we focused on binary NP structures, this work could easily be extended to ternary and quaternary NP superlattices. The additional NP species of distinct sizes or surface chemistry would introduce additional functionalities into the NP composite. However, fabrication of such multi-component NP superlattices is experimentally challenging and only a handful of ternary superlattices have so far been reported 19,43,59 . Recent MD simulations have shown that cubic NPs can be assembled into distinct 1D and 2D periodic architectures at fluid interfaces by taking advantage of their ability to tune the orientation of anisotropic NPs 60 . Superlattices of more complex-shaped NPs at interfaces could also be explored with our optimization strategy. The main challenge here would be the analytical treatment of interparticle and interfacial interactions, required for efficient discovery of assembly structures. Recent studies have indeed started to address these challenges, especially for nanocubes and other faceted NPs 61,62 . Our approach for searching for NP superlattices could be further improved by integrating strategies used in crystal structure prediction 63,64 , where lattice parameters, including angles and lengths of lattice vectors could be included as additional variables in the optimization. This may reduce the computational cost of optimization, as the large-area NP superlattices could be replaced by a single unit cell with periodic boundaries to avoid edge effects. While our model is most applicable to NPs much larger than the roughness of the interface and interparticle interactions dominated by vdW forces, our model can be readily extended to NP systems exhibiting other kinds of attractive interactions, such as those mediated by surface ligands or depletion forces. It is possible that the interactions between such NPs could deviate from the Lorentz-Berthelot rule assumed here, so the strength of inter-species attraction is no longer intermediate to the strengths of intra-species attraction, which could lead to other kinds of NP assemblies (Supplementary Fig. 17). In some cases, the surface of NPs is grafted with long polymer chains to assist in their assembly, where multibody effects become important 65 . Our optimization approach would then need to be combined with many-body potentials to uncover the more irregular, anisotropic NP assemblies expected to form as a result of multibody effects 21 .

Model for interparticle interactions
While NPs are typically grafted with molecules or polymer chains and interact with each other via a combination of solvent-dependent van der Waals (vdW) forces, solvent-mediated depletion forces, and graftmediated interactions, our model treats the two types of NPs as spheres of radii R 1 and R 2 that interact with each other, independent of the surrounding fluid, via a shifted Lennard-Jones (LJ) pairwise potential: where ij = 11, 22 and 12 denotes interactions between like and unlike pairs of NPs, r is the separation distance between the NPs, σ represents an intrinsic length scale of the system, and 4 ij = R i + R j À σ is a distance shift that prevents NPs from penetrating each other. We find that the resulting potential captures well the vdW interactions between small colloidal NPs ( Supplementary Fig. 1). The interaction strengths were assumed to be proportional to NP radii, consistent with the size scaling of vdW interactions between spherical colloids 38 , so that where we used Lorentz-Berthelot rules to infer the strength of interactions between unlike particles. The above model thus captures both the size-dependence in the depth of the interparticle interaction potential well and its narrower width (relative to particle size) compared to that of interatomic LJ potentials.

Model for interfacial free energy
We used an analytical model we recently developed 21 to describe the interfacial free energy of a surface-functionalized NP trapped at the interface ( Supplementary Fig. 2). According to the model, the interfacial free energies of the two species of NPs (relative to their reference bulk values in their favored fluid layer) as a function of their normal position z from the interface are given by where γ is the surface tension between the two fluids, and 4γ 1 γ 12 0 À γ 11 0 ≥ 0 and 4γ 2 γ 21 0 À γ 22 0 ≥ 0 are the differences in the surface energies of NP 1 and 2 with the two layers. When NPs intersect with interfacial plane (|z| < R i ), the free energy changes arise from two contributions: free energy gain due to the NP occluding part of the interface (first term) and free energy loss due to the unfavorable interactions between the NP and the incompatible layer (second term).
The equilibrium position z m of the NP (in the absence of interparticle interactions) can then be determined by ∂4F i =∂z| z m = 0, yielding Eq. (1) and (2).

Optimization approach
We used the basin-hopping Monte Carlo (BHMC) global optimization method for locating ground-state configurations of NPs at the fluidfluid interface [39][40][41] , that is, global minimum of the total energy of the system as given by Article https://doi.org/10.1038/s41467-022-35690-8 where r N represents the 3N coordinates of the NPs and I is a function that determines the identity of the particles (NP 1 or 2). In this algorithm, the original intricate potential energy surface E tot r N À Á is effectively transformed through local minimization into a staircase topography, represented by e E tot r N À Á = minfE tot r N À Á g, which is then explored via MC sampling using the Metropolis criterion. As shown in the flow chart (Fig. 1b), the initial seed configuration was randomly generated in a cuboid box in which the interfacial plane was fixed at z = 0. A sequence of MC moves was then implemented to produce new seed configurations, which were further optimized by using local minimization approaches. In particular, the Polak-Ribiere variant of the conjugate gradient algorithm 66 was used for local minimization and the inexact line search was performed based on the Wolfe conditions 67,68 . The Metropolis criterion was then used to update the NP configuration, which was then set as the seed structure for the next iteration. The above process was repeated until convergence was reached. While ground-state structures are strictly speaking only valid at zero temperature, differences in vibrational entropy of competing crystalline structures are often small. Thus, structures obtained via energy minimization may also represent minimum free-energy structures at non-zero temperatures.

New Monte Carlo moves
The MC displacement move implemented in standard BHMC algorithms for particle assembly is designed to search for 3D clusters and is not efficient for searching for the 2D structures typically obtained in interfacial assembly 13,19,21 . Inspired by the work of Rondina and Silva 41 , we introduced 5 additional MC moves into the standard BHMC algorithm for improving the sampling efficiency of NP configurations at interfaces (Fig. 1b): (i) Displacement move (D). This is the standard move employed in the original BHMC algorithm, which randomly displaces NPs around their original positions. The maximum displacements are fixed for all NPs. (ii) Position-based displacement move (P). In this variant of D, the magnitude of NP displacement depends on the distance of the NP from the centroid of the cluster 41,69 . This move takes advantage of differences in the energy landscape experienced by NPs at the surface versus core of the clusters. (iii) Twist move (T). This move separates NPs into two categories based on whether they are above or below the interfacial plane. The move then randomly selects one of the two categories and rotates all NPs in that category by a random angle about an axis normal to the interface passing through the cluster centroid. (iv) Interfacial move (I). It usually takes a long "diffusive time" for a NP fully immersed in one liquid to migrate toward the interface. Therefore, this move directly displaces the NP farthest from the interfacial plane to a plane near the interface. (v) Switch move (S). Since both species of NPs exhibit favorable interactions with one fluid layer and unfavorable interactions with the opposite layer, there is a free energy loss when NPs intrude into opposite layers. This move targets pairs of NPs of opposite species that intrude into their corresponding disfavored layer and then switches their type. (vi) Rotation move (R). This move chooses either the NP farthest from the cluster centroid or the NP with the least neighbors and randomly rotates the chosen NP on a circle parallel to the interfacial plane. This move is designed for filling in vacant sites in planar structures.
These moves were placed in a loop sequence of P-S-T-D-I-S-R-D-R. Each move in the sequence was executed for a period of MC steps until the move was consecutively rejected by the Metropolis criterion, upon which we switched to the next move in the sequence. We further realized that NPs could get trapped in very deep local minima from which they are unable to escape via the ongoing MC move, resulting in the same structure appearing again and again. To circumvent this, we implemented the ultrafast shape recognition algorithm 70,71 to identify duplicated NP structures. The BHMC algorithm would then directly jump to the next type of MC move if such duplicated structures are detected in consecutive MC steps. For optimization of large systems, an additional operation was added to the MC cycle, where the current best structure (till that MC step) was taken as the seed structure (instead of the instantaneous structure at that step) after every 2000 MC steps. This is because the current best structure was found to typically have more similarity to the ground-state configuration compared to the instantaneous structure.

Validation of optimization strategy
First, we adapted our optimization algorithm to the well-known bulk system of Lennard-Jones atom clusters, as there is no previous optimization work on 2D BNSLs at interfaces. Our algorithm was able to accurately locate the known global minima for three representative clusters containing up to 50 atoms 39,72,73 (Supplementary Fig. 18a). Second, we compared our predicted NP structures against observations from CGMD simulations that studied assembly of polymer-grafted NPs at a polymer-polymer interface 21 . These simulations showed that bare NPs that interact neutrally with both polymer layers assemble into a hexagonally arranged monolayer positioned symmetrically at the interface ( Supplementary Fig. 18b, bottom left); our optimization algorithm predicts the same morphology for a single species of neutral NPs, i.e., χ = 0 ( Supplementary Fig. 18b, top left). When co-assembling two species of grafted NPs with ligands that have different miscibility with the two fluids, a square bilayer was obtained in the CGMD simulations, which is reproduced by our algorithm for the case χ = 0:6 where the NPs strongly favor the opposite fluids of the interface (Supplementary Fig. 18b, middle). When the interactions between NPs (ϵ = 2:56) are large enough to overwhelm those between the NPs and the interface, the NPs form the same globular structure they form during bulk assembly; 74 this structure is also correctly reproduced by our optimization algorithm ( Supplementary Fig. 18b, right).

Phenomenological model of phase boundaries
The total free energy of a structural phase is composed of two terms: interfacial energy of the NPs and interparticle interaction energy. Using Eqs. (2) and (5), the interfacial energy (normalized by πR 2 1 γ) can be expressed as Àm 1 À χ ð Þ 2 , where m is the effective number of NPs at their equilibrium positions (z m ). As the interactions between NPs are short ranged, the interaction energy between NPs (normalized by πR 2 1 γ) can be approximated by Ànϵ (see Eq. 4), where n is the total number of NP-NP contacts exhibited by the structure. Therefore, the total free energy of a structure (normalized by πR 2 1 γ) can be represented by Àm 1 À χ ð Þ 2 À nϵ . For monolayer structures in the phase diagram, all particles can sit at their preferred interfacial positions without sacrificing any interparticle interactions, so m = N = 12, and counting NP-NP contacts in the structure reveals that n = 24 ( Fig. 3a; Supplementary Table 1). For the bilayer structures, m can be determined from the ratio of their total interfacial energy and the minimum value of the interfacial free energy, that is, À 1 À χ ð Þ 2 . Based on geometry, the half-spacing between the two layers of NPs in the square and hexagonal bilayers is given by z s = 2:16σ and 2:5σ. These values can be substituted into the Eqs. 9 and 10 to obtain the interfacial energy, and thereby m. We find that m≈11 for the square bilayer and m≈10 for the hexagonal bilayer at χ = 0:6. Furthermore, we obtain n = 30 for the square bilayer and n = 31 for the hexagonal bilayer ( Fig. 3b; Supplementary Table 1). Lastly, for the globular clusters, specifically the truncated triangular dipyramid structure, n = 33. Since the globular clusters can adjust their position and orientation with respect to the interface plane as χ varies, m could change with χ. To this end, we take an averaged value of m for all structures above or below χ = 0:5, which give m ≈ 6:5 for χ > 0:5 and m ≈ 4:5 for χ<0:5. (Fig. 3c; Supplementary Table 1). By equating the total energies from pairs of the structural phases, we obtain the following three equations for boundaries between the three structural phases (Supplementary Fig. 5): • Boundary between monolayers and bilayers (blue line): • Boundary between bilayers and globular clusters (red line): ϵ = À2:18 1 À χ ð Þ 2 À 12 1 À 0:833 ð Þ χ + 6ð1 À 0:833 2 Þ ð13Þ • Boundary between monolayers and globular clusters (purple line): Predicting transition between triangular and rhombic monolayers Since the interaction strength between NPs is proportional to particle size, the unit cell of a triangular lattice (with effectively one contact formed between large NPs) would display more favorable interparticle interaction energy than that of a rhombic lattice (with one contact formed between two small NP). According to Eq. (8), this interaction energy difference is roughly equal to the difference ε 11 À ε 22 = ε 11 ð1 À R 2 =R 1 Þ, yielding a dimensionless energy of 2=3ϵ via Eq. (4). On the other hand, the triangular monolayer displays less favorable interfacial energy than its rhombic counterpart, as the small NPs exhibit more deviation from their ideal position ( z h i≈0:8σ compared to z m = χR 2 = 0:4σ) in the triangular monolayer ( Supplementary  Fig. 14d) than in the rhombic monolayer, where the small NPs are located almost at the ideal position ( Supplementary Fig. 14e). In both monolayers, the large NPs reside very close to their ideal position (z m = χR 1 = 1:2σ), as large NPs sacrifice much larger interfacial energy compared to the small ones due to their × 9 larger cross-sectional area. The interfacial energy difference between the unit cells of the two monolayers can then be estimated by the difference in the interfacial energy of their two small NPs, as given by 2½4F 2 0:8σ ð ÞÀ4F 2 ð0:4σÞ, which yields a dimensionless energy value of 0:0356 via Eq. 10. When ϵ > 0:053 (obtained by setting 2=3ϵ > 0:0356), the loss in interfacial energy can be compensated by the gain in interparticle interactions, favoring the formation of the triangular monolayer (Fig. 4d). Otherwise, the NPs prefer to maximize their interfacial energy by forming the rhombic monolayer (Fig. 4e).

Model for predicting vdW interaction energy of spherical NPs
Assuming that the atoms in the two NPs are of size σ and interact via the LJ potential, the vdW energy W ðDÞ between two NPs of radii R separated by a surface-to-surface distance of D ð≪RÞ is given by W D ð Þ= ARðσ 6 =2520D 7 À 1=12DÞ 38 . At the energy minimum (30 À1=6 σ), the vdW attraction is given by ≈ À0:126AR=σ. Since the Hamaker constant of silica NPs (σ = 0.21 nm) in air is 65 zJ 75 , and that of silver NPs (σ = 0.172 nm) in water is 61.6 zJ 58 , the interparticle interaction strength between silica NPs and between silver NPs are given by ε 22 = 117 zJ and ε 11 = 135 zJ, using R = 3 nm.

Data availability
The data generated in this study have been deposited in the public GitHub (https://github.com/yilongz393/Tailored-BHMC) without any restrictions. Source data are provided with this paper.

Code availability
The codes used in this study have been deposited in the public GitHub (https://github.com/yilongz393/ Tailored-BHMC) without any restrictions.